####rHAPC Models (Secondary model, interviewed only)######

model_CCREM_Fullcountry_onlyInterviewed = datCompleteNew_interviewed %>% 
  filter(waves %in% c("EVS1","EVS2","EVS3","EVS4","EVS5")) %>% 
  filter(!(dv_name == "democracy_freeElection" | dv_name == "democracy_ArmyTakeOver" | dv_name == "democracy_protectLiberty" | dv_name == "DemImportant")) %>% 
  group_by(dv_name, country) %>%
  do(fitCountry = lmer(dv_value ~ age_grp_centered + I(age_grp_centered^2) + I(age_grp_centered^3) + born_adult_centered  + I(born_adult_centered^2)+ I(born_adult_centered^3) + sex + (1|waves.f) + (1|age_cohort_five), data = .)
  )

#Save the regression output for the Shiny App
source("16_RegTables_rHAPC_onlyInterviewed.R")

name_flag <- model_CCREM_Fullcountry_onlyInterviewed$dv_name
country_flag <- model_CCREM_Fullcountry_onlyInterviewed$country


#Create new data for predictions
datBind <- list()

for(i in 1:length(name_flag)){
  datBind[[i]] <- datCompleteNew_interviewed %>% 
    filter(waves %in% c("EVS1","EVS2","EVS3","EVS4","EVS5")) %>% 
    filter(!(dv_name == "democracy_freeElection" | dv_name == "democracy_ArmyTakeOver" | dv_name == "democracy_protectLiberty" | dv_name == "DemImportant")) %>% 
    filter(dv_name == name_flag[i]) %>% 
    filter(country == country_flag[i]) %>% 
    drop_na(age_grp_centered, waves.f, born_adult_centered, age_cohort_five, dv_value, sex) %>% 
    dplyr::select(dv_name, dv_value, age_grp_centered, waves.f, born_adult_centered, born_adult, age_grp, age_cohort_five, sex) %>% 
    droplevels() %>% 
    dplyr::select(age_grp_centered, waves.f, born_adult_centered, born_adult, age_grp, age_cohort_five, sex)
}

predList <- list()

for(i in 1:length(name_flag)){
  predList[[i]] <- predict(model_CCREM_Fullcountry_onlyInterviewed$fitCountry[[i]], datBind[[i]],  re.form = NULL)
}

mmList <- list()

for(i in 1:length(name_flag)){
  mmList[[i]] <- model.matrix(~age_grp_centered + I(age_grp_centered^2) + I(age_grp_centered^3) + born_adult_centered + I(born_adult_centered^2) + I(born_adult_centered^3) + sex,datBind[[i]])
}

predFun<-function(.) mmList[[i]]%*%fixef(.)

bbList <- list()
bb_seList <- list()

for(i in 1:length(name_flag)){
  bbList[[i]] <- bootMer(model_CCREM_Fullcountry_onlyInterviewed$fitCountry[[i]],FUN=predFun,nsim=200, re.form = NULL)
  bb_seList[[i]] <- apply(bbList[[i]]$t,2,function(x) x[order(x)][c(5,195)])
  print(i)
}

for(i in 1:length(name_flag)){
  datBind[[i]]$LL <- bb_seList[[i]][1,]
  datBind[[i]]$UL <- bb_seList[[i]][2,]
}

preds <- list()

for(i in 1:length(name_flag)){
  preds[[i]] <- transform(cbind(data.frame(predList[[i]]), datBind[[i]]))
}




#Aggregate predictions
predBornAdult <- list()

for(i in 1:length(name_flag)){
  
  predBornAdult[[i]] <- preds[[i]] %>% 
    mutate(se.fit = (`predList..i..` - LL)/1.96,
           fit = `predList..i..`) %>% 
    group_by(born_adult) %>% 
    mutate(m = n(),
           mean_fit = mean(fit),
           S_sq = (fit - mean_fit)^2/(m - 1)) %>%
    replace_na(list(S_sq = 0)) %>% 
    summarize(fit = mean(fit),
              S_sq = mean(S_sq),
              se.fit = mean(se.fit^2) + S_sq,
              se.fit = sqrt(se.fit)) %>% 
    ungroup() %>% 
    mutate(UL = fit + 1.96 * se.fit,
           LL = fit - 1.96 * se.fit) %>% 
    dplyr::select(born_adult, fit, UL, LL)
  
}



predAgeGrp <- list()

for(i in 1:length(name_flag)){
  
  predAgeGrp[[i]] <- preds[[i]] %>% 
    mutate(se.fit = (`predList..i..` - LL)/1.96,
           fit = `predList..i..`) %>% 
    group_by(age_grp) %>% 
    mutate(m = n(),
           mean_fit = mean(fit),
           S_sq = (fit - mean_fit)^2/(m - 1)) %>%
    replace_na(list(S_sq = 0)) %>% 
    summarize(fit = mean(fit),
              S_sq = mean(S_sq),
              se.fit = mean(se.fit^2) + S_sq,
              se.fit = sqrt(se.fit)) %>% 
    ungroup() %>% 
    mutate(UL = fit + 1.96 * se.fit,
           LL = fit - 1.96 * se.fit) %>% 
    dplyr::select(age_grp, fit, UL, LL)
  
}




predPeriod <- list()

for(i in 1:length(name_flag)){
  
  predPeriod[[i]] <- preds[[i]] %>% 
    mutate(se.fit = (`predList..i..` - LL)/1.96,
           fit = `predList..i..`) %>% 
    group_by(waves.f) %>% 
    mutate(m = n(),
           mean_fit = mean(fit),
           S_sq = (fit - mean_fit)^2/(m - 1)) %>%
    replace_na(list(S_sq = 0)) %>% 
    summarize(fit = mean(fit),
              S_sq = mean(S_sq),
              se.fit = mean(se.fit^2) + S_sq,
              se.fit = sqrt(se.fit)) %>% 
    ungroup() %>% 
    mutate(UL = fit + 1.96 * se.fit,
           LL = fit - 1.96 * se.fit) %>% 
    dplyr::select(waves.f, fit, UL, LL)
  
}


names(predBornAdult) <- c(paste(name_flag, as.character(country_flag), sep = ";"))

predBornAdult <- bind_rows(predBornAdult, .id = "dv_name")

predBornAdult <- predBornAdult %>% 
  tidyr::separate(dv_name, c("dv_name", "country"), ";")

names(predAgeGrp) <- c(paste(name_flag, as.character(country_flag), sep = ";"))

predAgeGrp <- bind_rows(predAgeGrp, .id = "dv_name")

predAgeGrp <- predAgeGrp %>% 
  tidyr::separate(dv_name, c("dv_name", "country"), ";")


names(predPeriod) <- c(paste(name_flag, as.character(country_flag), sep = ";"))

predPeriod <- bind_rows(predPeriod, .id = "dv_name")

predPeriod <- predPeriod %>% 
  tidyr::separate(dv_name, c("dv_name", "country"), ";") %>% 
  mutate(year = recode(waves.f, "EVS1" = 1981, "EVS2" = 1990, "EVS3" = 1999, "EVS4" = 2008, "EVS5" = 2017))





#Create plot of predictions for all countries and all DV and store them as png
pg_list <- list()

for(i in 1:length(predList)){
  
  message(i)
  
  pg_list[[i]] <- local({
    i <- i 
    pp.data <- predBornAdult %>% filter(dv_name == name_flag[i], country == country_flag[i])
    pp_period.data <- predPeriod %>% filter(dv_name == name_flag[i], country == country_flag[i])
    p1 <- ggplot() + 
      geom_errorbar(data = pp.data, aes(x=born_adult, ymin = LL, ymax = UL), alpha = 0.2, colour = cols[2]) +
      geom_point(data = pp.data, aes(x=born_adult, y=fit), alpha = 0.1, size = 1, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = fit), se = F, alpha = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = LL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = UL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) +
      ggtitle(paste(pp.data$country[1])) + 
      scale_y_continuous(name="", breaks=seq(-1, 1, 0.4), limits = c(-1,1)) + 
      geom_point(data = pp_period.data, aes(x = (((year - 1980)*100)/40) + 1920, y = fit), color = cols[1]) + 
      geom_text(data = pp_period.data, aes(x = (((year - 1980)*100)/40) + 1920, y = fit, label=year),vjust = 1, hjust = 1, size = 2) + 
      scale_x_continuous(name = " ",breaks = c(1920, 1940, 1960, 1980, 2000, 2020),limits = c(1920, 2020), sec.axis = sec_axis(~(((. - 1920)*40)/100) + 1980, breaks = c(1980, 2000, 2020), name = "") ) + 
      theme(
        axis.line = element_line(colour = "black"), 
        text = element_text(size=7), 
        axis.text = element_text(size = 7),
        plot.title = element_text(margin = margin(0,0,0,0), size = 7, vjust = -4),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt")
      ) + 
      theme_light() + 
      theme(
        plot.title = element_text(margin = margin(0,0,-5,0)),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt"),
        axis.text.x.bottom = element_text(angle = -90, hjust = 0, vjust = 0)
      )
    print(p1)
  })
  
  
}

for(i in 1:24){
  pg_list[[i]] <- pg_list[[i]] + scale_y_continuous(name="", limits = c(-1.1,1.3)) 
}


for(i in 25:length(predList)){
  pg_list[[i]] <- pg_list[[i]] + scale_y_continuous(name="", breaks=seq(0, 1, 0.2), limits = c(-0.1,1.1)) 
}


for(i in seq(1,length(predList), by = 24)){
  agg_png(file=paste0(wdShinyPlots, "plotsFull_CCREM_cohort_",i,"_onlyInterviewed.png"), width = 3360, height = 2240, res = 300)
  
  grid.arrange(arrangeGrob(pg_list[[i]], pg_list[[i+1]], pg_list[[i+2]], pg_list[[i+3]], pg_list[[i+4]], pg_list[[i+5]], pg_list[[i+6]], pg_list[[i+7]], pg_list[[i+8]], pg_list[[i+9]], pg_list[[i+10]], pg_list[[i+11]], pg_list[[i+12]], pg_list[[i+13]], pg_list[[i+14]], pg_list[[i+15]], pg_list[[i+16]], pg_list[[i+17]], pg_list[[i+18]], pg_list[[i+19]], pg_list[[i+20]], pg_list[[i+21]], pg_list[[i+22]], pg_list[[i+23]], nrow = 4))
  
  dev.off()
}





plc_list <- list()

for(i in 1:length(predList)){
  
  message(i)
  
  plc_list[[i]] <- local({
    i <- i 
    pp.data <- predAgeGrp %>% filter(dv_name == name_flag[i], country == country_flag[i])
    p1 <- ggplot() + 
      geom_errorbar(data = pp.data, aes(x=age_grp, ymin = LL, ymax = UL), alpha = 0.2, colour = cols[2]) +
      geom_point(data = pp.data, aes(x=age_grp, y=fit), alpha = 0.1, size = 1, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=age_grp,y = fit), se = F, alpha = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=age_grp,y = LL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=age_grp,y = UL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) +
      ggtitle(paste(pp.data$country[1])) + 
      scale_y_continuous(name="", breaks=seq(-1, 1, 0.4), limits = c(-1,1)) + 
      scale_x_continuous(name = " ",breaks = seq(1,8,1) ,labels = c("15-19", "", "30-39", "", "50-59", "", "70-79", "")) + 
      theme(
        axis.line = element_line(colour = "black"), 
        text = element_text(size=7), 
        axis.text = element_text(size = 7),
        plot.title = element_text(margin = margin(0,0,0,0), size = 7, vjust = -4),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt")
      ) + 
      theme_light() + 
      theme(
        plot.title = element_text(margin = margin(0,0,0,0)),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt"),
        axis.text.x.bottom = element_text(angle = -90, hjust = 0, vjust = 0)
      )
    print(p1)
  })
  
  
}

for(i in 25:length(predList)){
  plc_list[[i]] <- plc_list[[i]] + scale_y_continuous(name="", breaks=seq(0, 1, 0.2), limits = c(-0.1,1.1)) 
}


for(i in seq(1,length(predList), by = 24)){
  agg_png(file=paste0(wdShinyPlots, "plotsFull_CCREM_lifecycle_",i,"_onlyInterviewed.png"), width = 3360, height = 2240, res = 300)
  
  grid.arrange(arrangeGrob(plc_list[[i]], plc_list[[i+1]], plc_list[[i+2]], plc_list[[i+3]], plc_list[[i+4]], plc_list[[i+5]], plc_list[[i+6]], plc_list[[i+7]], plc_list[[i+8]], plc_list[[i+9]], plc_list[[i+10]], plc_list[[i+11]], plc_list[[i+12]], plc_list[[i+13]], plc_list[[i+14]], plc_list[[i+15]], plc_list[[i+16]], plc_list[[i+17]], plc_list[[i+18]], plc_list[[i+19]], plc_list[[i+20]], plc_list[[i+21]], plc_list[[i+22]], plc_list[[i+23]], nrow = 4))
  
  dev.off()
}




#Run the GAM models for the smaller country sample (Democracy Important & Essential variables)
model_CCREM_Subcountry_DemImportant_onlyInterviewed = datCompleteNew_interviewed %>% 
  filter(country %in% countrySub_DemImportant) %>% 
  filter((dv_name == "DemImportant")) %>% 
  group_by(dv_name, country) %>%
  do(fitCountry = lmer(dv_value ~ age_grp_centered + I(age_grp_centered^2) + I(age_grp_centered^3) + born_adult_centered  + I(born_adult_centered^2)+ I(born_adult_centered^3) + sex + (1|waves.f) + (1|age_cohort_five), data = .)
  )

model_CCREM_Subcountry_Essentials_onlyInterviewed = datCompleteNew_interviewed %>% 
  filter(country %in% countrySub_Essentials) %>% 
  filter((dv_name == "democracy_freeElection" | dv_name == "democracy_ArmyTakeOver" | dv_name == "democracy_protectLiberty")) %>% 
  group_by(dv_name, country) %>%
  do(fitCountry = lmer(dv_value ~ age_grp_centered + I(age_grp_centered^2) + I(age_grp_centered^3) + born_adult_centered  + I(born_adult_centered^2)+ I(born_adult_centered^3) + sex + (1|waves.f) + (1|age_cohort_five), data = .)
  )

model_CCREM_Subcountry_onlyInterviewed <- bind_rows(model_CCREM_Subcountry_DemImportant_onlyInterviewed, model_CCREM_Subcountry_Essentials_onlyInterviewed)

#Save the regression output for the Shiny App
source("18_RegTables_rHAPC_Subsample_onlyInterviewed.R")



name_flag_DemImportant <- model_CCREM_Subcountry_DemImportant_onlyInterviewed$dv_name
name_flag_Essentials <- model_CCREM_Subcountry_Essentials_onlyInterviewed$dv_name

country_flag_DemImportant <- model_CCREM_Subcountry_DemImportant_onlyInterviewed$country
country_flag_Essentials <- model_CCREM_Subcountry_Essentials_onlyInterviewed$country


#Create new data for predictions
datBind_DemImportant <- list()
datBind_Essentials <- list()



for(i in 1:length(name_flag_DemImportant)){
  datBind_DemImportant[[i]] <- datCompleteNew_interviewed %>% 
    filter(waves %in% c("EVS5", "WVS5", "WVS6")) %>% 
    filter(dv_name == "DemImportant") %>% 
    filter(dv_name == name_flag_DemImportant[i]) %>% 
    filter(country == country_flag_DemImportant[i]) %>% 
    drop_na(age_grp_centered, waves.f, born_adult_centered, age_cohort_five, dv_value, sex) %>% 
    dplyr::select(dv_name, dv_value, age_grp_centered, waves.f, born_adult_centered, born_adult, age_grp, age_cohort_five, sex) %>% 
    droplevels() %>% 
    dplyr::select(age_grp_centered, waves.f, born_adult_centered, born_adult, age_grp, age_cohort_five, sex)
}

for(i in 1:length(name_flag_Essentials)){
  datBind_Essentials[[i]] <- datCompleteNew_interviewed %>% 
    filter(waves %in% c("EVS5", "WVS5", "WVS6")) %>% 
    filter(dv_name == "democracy_freeElection" | dv_name == "democracy_ArmyTakeOver" | dv_name == "democracy_protectLiberty") %>% 
    filter(dv_name == name_flag_Essentials[i]) %>% 
    filter(country == country_flag_Essentials[i]) %>% 
    drop_na(age_grp_centered, waves.f, born_adult_centered, age_cohort_five, dv_value, sex) %>% 
    dplyr::select(dv_name, dv_value, age_grp_centered, waves.f, born_adult_centered, born_adult, age_grp, age_cohort_five, sex) %>% 
    droplevels() %>% 
    dplyr::select(age_grp_centered, waves.f, born_adult_centered, born_adult, age_grp, age_cohort_five, sex)
  
}





predList_DemImportant <- list()
predList_Essentials <- list()

for(i in 1:length(name_flag_DemImportant)){
  predList_DemImportant[[i]] <- predict(model_CCREM_Subcountry_DemImportant_onlyInterviewed$fitCountry[[i]], datBind_DemImportant[[i]],  re.form = NULL)
}

for(i in 1:length(name_flag_Essentials)){
  predList_Essentials[[i]] <- predict(model_CCREM_Subcountry_Essentials_onlyInterviewed$fitCountry[[i]], datBind_Essentials[[i]],  re.form = NULL)
}

mmList_DemImportant <- list()
mmList_Essentials <- list()

for(i in 1:length(name_flag_DemImportant)){
  mmList_DemImportant[[i]] <- model.matrix(~age_grp_centered + I(age_grp_centered^2) + I(age_grp_centered^3) + born_adult_centered + I(born_adult_centered^2) + I(born_adult_centered^3) + sex, datBind_DemImportant[[i]])
}

for(i in 1:length(name_flag_Essentials)){
  mmList_Essentials[[i]] <- model.matrix(~age_grp_centered + I(age_grp_centered^2) + I(age_grp_centered^3) + born_adult_centered + I(born_adult_centered^2) + I(born_adult_centered^3) + sex, datBind_Essentials[[i]])
}

predFun_DemImportant <- function(.) mmList_DemImportant[[i]]%*%fixef(.)
predFun_Essentials <- function(.) mmList_Essentials[[i]]%*%fixef(.)

bbList_DemImportant <- list()
bb_seList_DemImportant <- list()

bbList_Essentials <- list()
bb_seList_Essentials <- list()

for(i in 1:length(name_flag_DemImportant)){
  bbList_DemImportant[[i]] <- bootMer(model_CCREM_Subcountry_DemImportant_onlyInterviewed$fitCountry[[i]],FUN=predFun_DemImportant,nsim=200, re.form = NULL)
  bb_seList_DemImportant[[i]] <- apply(bbList_DemImportant[[i]]$t,2,function(x) x[order(x)][c(5,195)])
  print(i)
}

for(i in 1:length(name_flag_Essentials)){
  bbList_Essentials[[i]] <- bootMer(model_CCREM_Subcountry_Essentials_onlyInterviewed$fitCountry[[i]],FUN=predFun_Essentials,nsim=200, re.form = NULL)
  bb_seList_Essentials[[i]] <- apply(bbList_Essentials[[i]]$t,2,function(x) x[order(x)][c(5,195)])
  print(i)
}

for(i in 1:length(name_flag_DemImportant)){
  datBind_DemImportant[[i]]$LL <- bb_seList_DemImportant[[i]][1,]
  datBind_DemImportant[[i]]$UL <- bb_seList_DemImportant[[i]][2,]
}

for(i in 1:length(name_flag_Essentials)){
  datBind_Essentials[[i]]$LL <- bb_seList_Essentials[[i]][1,]
  datBind_Essentials[[i]]$UL <- bb_seList_Essentials[[i]][2,]
}


preds_DemImportant <- list()
preds_Essentials <- list()

for(i in 1:length(name_flag_DemImportant)){
  preds_DemImportant[[i]] <- transform(cbind(data.frame(predList_DemImportant[[i]]), datBind_DemImportant[[i]]))
}

for(i in 1:length(name_flag_Essentials)){
  preds_Essentials[[i]] <- transform(cbind(data.frame(predList_Essentials[[i]]), datBind_Essentials[[i]]))
}




#Aggregate predictions
predBornAdult_DemImportant <- list()
predBornAdult_Essentials <- list()

for(i in 1:length(name_flag_DemImportant)){
  
  predBornAdult_DemImportant[[i]] <- preds_DemImportant[[i]] %>% 
    mutate(se.fit = (`predList_DemImportant..i..` - LL)/1.96,
           fit = `predList_DemImportant..i..`) %>% 
    group_by(born_adult) %>% 
    mutate(m = n(),
           mean_fit = mean(fit),
           S_sq = (fit - mean_fit)^2/(m - 1)) %>%
    replace_na(list(S_sq = 0)) %>% 
    summarize(fit = mean(fit),
              S_sq = mean(S_sq),
              se.fit = mean(se.fit^2) + S_sq,
              se.fit = sqrt(se.fit)) %>% 
    ungroup() %>% 
    mutate(UL = fit + 1.96 * se.fit,
           LL = fit - 1.96 * se.fit) %>% 
    dplyr::select(born_adult, fit, UL, LL)
  
}

for(i in 1:length(name_flag_Essentials)){
  
  predBornAdult_Essentials[[i]] <- preds_Essentials[[i]] %>% 
    mutate(se.fit = (`predList_Essentials..i..` - LL)/1.96,
           fit = `predList_Essentials..i..`) %>% 
    group_by(born_adult) %>% 
    mutate(m = n(),
           mean_fit = mean(fit),
           S_sq = (fit - mean_fit)^2/(m - 1)) %>%
    replace_na(list(S_sq = 0)) %>% 
    summarize(fit = mean(fit),
              S_sq = mean(S_sq),
              se.fit = mean(se.fit^2) + S_sq,
              se.fit = sqrt(se.fit)) %>% 
    ungroup() %>% 
    mutate(UL = fit + 1.96 * se.fit,
           LL = fit - 1.96 * se.fit) %>% 
    dplyr::select(born_adult, fit, UL, LL)
  
}



predAgeGrp_DemImportant <- list()
predAgeGrp_Essentials <- list()

for(i in 1:length(name_flag_DemImportant)){
  
  predAgeGrp_DemImportant[[i]] <- preds_DemImportant[[i]] %>% 
    mutate(se.fit = (`predList_DemImportant..i..` - LL)/1.96,
           fit = `predList_DemImportant..i..`) %>% 
    group_by(age_grp) %>% 
    mutate(m = n(),
           mean_fit = mean(fit),
           S_sq = (fit - mean_fit)^2/(m - 1)) %>%
    replace_na(list(S_sq = 0)) %>% 
    summarize(fit = mean(fit),
              S_sq = mean(S_sq),
              se.fit = mean(se.fit^2) + S_sq,
              se.fit = sqrt(se.fit)) %>% 
    ungroup() %>% 
    mutate(UL = fit + 1.96 * se.fit,
           LL = fit - 1.96 * se.fit) %>% 
    dplyr::select(age_grp, fit, UL, LL)
  
}


for(i in 1:length(name_flag_Essentials)){
  
  predAgeGrp_Essentials[[i]] <- preds_Essentials[[i]] %>% 
    mutate(se.fit = (`predList_Essentials..i..` - LL)/1.96,
           fit = `predList_Essentials..i..`) %>% 
    group_by(age_grp) %>% 
    mutate(m = n(),
           mean_fit = mean(fit),
           S_sq = (fit - mean_fit)^2/(m - 1)) %>%
    replace_na(list(S_sq = 0)) %>% 
    summarize(fit = mean(fit),
              S_sq = mean(S_sq),
              se.fit = mean(se.fit^2) + S_sq,
              se.fit = sqrt(se.fit)) %>% 
    ungroup() %>% 
    mutate(UL = fit + 1.96 * se.fit,
           LL = fit - 1.96 * se.fit) %>% 
    dplyr::select(age_grp, fit, UL, LL)
  
}




predPeriod_DemImportant <- list()
predPeriod_Essentials <- list()

for(i in 1:length(name_flag_DemImportant)){
  
  predPeriod_DemImportant[[i]] <- preds_DemImportant[[i]] %>% 
    mutate(se.fit = (`predList_DemImportant..i..` - LL)/1.96,
           fit = `predList_DemImportant..i..`) %>% 
    group_by(waves.f) %>% 
    mutate(m = n(),
           mean_fit = mean(fit),
           S_sq = (fit - mean_fit)^2/(m - 1)) %>%
    replace_na(list(S_sq = 0)) %>% 
    summarize(fit = mean(fit),
              S_sq = mean(S_sq),
              se.fit = mean(se.fit^2) + S_sq,
              se.fit = sqrt(se.fit)) %>% 
    ungroup() %>% 
    mutate(UL = fit + 1.96 * se.fit,
           LL = fit - 1.96 * se.fit) %>% 
    dplyr::select(waves.f, fit, UL, LL)
  
}

for(i in 1:length(name_flag_Essentials)){
  
  predPeriod_Essentials[[i]] <- preds_Essentials[[i]] %>% 
    mutate(se.fit = (`predList_Essentials..i..` - LL)/1.96,
           fit = `predList_Essentials..i..`) %>% 
    group_by(waves.f) %>% 
    mutate(m = n(),
           mean_fit = mean(fit),
           S_sq = (fit - mean_fit)^2/(m - 1)) %>%
    replace_na(list(S_sq = 0)) %>% 
    summarize(fit = mean(fit),
              S_sq = mean(S_sq),
              se.fit = mean(se.fit^2) + S_sq,
              se.fit = sqrt(se.fit)) %>% 
    ungroup() %>% 
    mutate(UL = fit + 1.96 * se.fit,
           LL = fit - 1.96 * se.fit) %>% 
    dplyr::select(waves.f, fit, UL, LL)
  
}


names(predBornAdult_DemImportant) <- c(paste(name_flag_DemImportant, as.character(country_flag_DemImportant), sep = ";"))
names(predBornAdult_Essentials) <- c(paste(name_flag_Essentials, as.character(country_flag_Essentials), sep = ";"))

predBornAdult_DemImportant <- bind_rows(predBornAdult_DemImportant, .id = "dv_name")
predBornAdult_Essentials <- bind_rows(predBornAdult_Essentials, .id = "dv_name")

predBornAdult_DemImportant <- predBornAdult_DemImportant %>% 
  tidyr::separate(dv_name, c("dv_name", "country"), ";")
predBornAdult_Essentials <- predBornAdult_Essentials %>% 
  tidyr::separate(dv_name, c("dv_name", "country"), ";")

names(predAgeGrp_DemImportant) <- c(paste(name_flag_DemImportant, as.character(country_flag_DemImportant), sep = ";"))
names(predAgeGrp_Essentials) <- c(paste(name_flag_Essentials, as.character(country_flag_Essentials), sep = ";"))

predAgeGrp_DemImportant <- bind_rows(predAgeGrp_DemImportant, .id = "dv_name")
predAgeGrp_Essentials <- bind_rows(predAgeGrp_Essentials, .id = "dv_name")

predAgeGrp_DemImportant <- predAgeGrp_DemImportant %>% 
  tidyr::separate(dv_name, c("dv_name", "country"), ";") %>% 
  mutate(age_grp2 = recode(age_grp, "1" = "15-19", "2" = "20-29", "3" = "30-39", "4" = "40-49", "5" = "50-59", "6" = "60-69", "7" = "70-79", "8" = "80+" ))
predAgeGrp_Essentials <- predAgeGrp_Essentials %>% 
  tidyr::separate(dv_name, c("dv_name", "country"), ";") %>% 
  mutate(age_grp2 = recode(age_grp, "1" = "15-19", "2" = "20-29", "3" = "30-39", "4" = "40-49", "5" = "50-59", "6" = "60-69", "7" = "70-79", "8" = "80+" ))


names(predPeriod_DemImportant) <- c(paste(name_flag_DemImportant, as.character(country_flag_DemImportant), sep = ";"))
names(predPeriod_Essentials) <- c(paste(name_flag_Essentials, as.character(country_flag_Essentials), sep = ";"))

predPeriod_DemImportant <- bind_rows(predPeriod_DemImportant, .id = "dv_name")
predPeriod_Essentials <- bind_rows(predPeriod_Essentials, .id = "dv_name")

predPeriod_DemImportant <- predPeriod_DemImportant %>% 
  tidyr::separate(dv_name, c("dv_name", "country"), ";") %>% 
  mutate(year = recode(waves.f, "WVS5" = 2007, "WVS6" = 2012, "EVS5" = 2017))
predPeriod_Essentials <- predPeriod_Essentials %>% 
  tidyr::separate(dv_name, c("dv_name", "country"), ";") %>% 
  mutate(year = recode(waves.f, "WVS5" = 2007, "WVS6" = 2012, "EVS5" = 2017))

#####Plotting######

pg_list_DemImportant <- list()
pg_list_Essentials <- list()


for(i in 1:length(predList_DemImportant)){
  
  message(i)
  
  pg_list_DemImportant[[i]] <- local({
    i <- i 
    pp.data <- predBornAdult_DemImportant %>% filter(dv_name == name_flag_DemImportant[i], country == country_flag_DemImportant[i])
    pp_period.data <- predPeriod_DemImportant %>% filter(dv_name == name_flag_DemImportant[i], country == country_flag_DemImportant[i])
    p1 <- ggplot() + 
      geom_errorbar(data = pp.data, aes(x=born_adult, ymin = LL, ymax = UL), alpha = 0.2, colour = cols[2]) +
      geom_point(data = pp.data, aes(x=born_adult, y=fit), alpha = 0.1, size = 1, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = fit), se = F, alpha = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = LL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = UL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) +
      ggtitle(paste(pp.data$country[1])) + 
      scale_y_continuous(name="", breaks=seq(-1, 1, 0.4), limits = c(-1,1)) + 
      geom_point(data = pp_period.data, aes(x = (((year - 1980)*100)/40) + 1920, y = fit), color = cols[1]) + 
      geom_text(data = pp_period.data, aes(x = (((year - 1980)*100)/40) + 1920, y = fit, label=year),vjust = 1, hjust = 1, size = 2) + 
      scale_x_continuous(name = " ",breaks = c(1920, 1940, 1960, 1980, 2000, 2020),limits = c(1920, 2020), sec.axis = sec_axis(~(((. - 1920)*40)/100) + 1980, breaks = c(1980, 2000, 2020), name = "") ) + 
      theme(
        axis.line = element_line(colour = "black"), 
        text = element_text(size=7), 
        axis.text = element_text(size = 7),
        plot.title = element_text(margin = margin(0,0,0,0), size = 7, vjust = -4),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt")
      ) + 
      theme_light() + 
      theme(
        plot.title = element_text(margin = margin(0,0,-5,0)),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt"),
        axis.text.x.bottom = element_text(angle = -90, hjust = 0, vjust = 0)
      )
    print(p1)
  })
  
  
}


for(i in 1:length(predList_Essentials)){
  
  message(i)
  
  pg_list_Essentials[[i]] <- local({
    i <- i 
    pp.data <- predBornAdult_Essentials %>% filter(dv_name == name_flag_Essentials[i], country == country_flag_Essentials[i])
    pp_period.data <- predPeriod_Essentials %>% filter(dv_name == name_flag_Essentials[i], country == country_flag_Essentials[i])
    p1 <- ggplot() + 
      geom_errorbar(data = pp.data, aes(x=born_adult, ymin = LL, ymax = UL), alpha = 0.2, colour = cols[2]) +
      geom_point(data = pp.data, aes(x=born_adult, y=fit), alpha = 0.1, size = 1, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = fit), se = F, alpha = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = LL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = UL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) +
      ggtitle(paste(pp.data$country[1])) + 
      scale_y_continuous(name="", breaks=seq(-1, 1, 0.4), limits = c(-1,1)) + 
      geom_point(data = pp_period.data, aes(x = (((year - 1980)*100)/40) + 1920, y = fit), color = cols[1]) + 
      geom_text(data = pp_period.data, aes(x = (((year - 1980)*100)/40) + 1920, y = fit, label=year),vjust = 1, hjust = 1, size = 2) + 
      scale_x_continuous(name = " ",breaks = c(1920, 1940, 1960, 1980, 2000, 2020),limits = c(1920, 2020), sec.axis = sec_axis(~(((. - 1920)*40)/100) + 1980, breaks = c(1980, 2000, 2020), name = "") ) + 
      theme(
        axis.line = element_line(colour = "black"), 
        text = element_text(size=7), 
        axis.text = element_text(size = 7),
        plot.title = element_text(margin = margin(0,0,0,0), size = 7, vjust = -4),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt")
      ) + 
      theme_light() + 
      theme(
        plot.title = element_text(margin = margin(0,0,-5,0)),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt"),
        axis.text.x.bottom = element_text(angle = -90, hjust = 0, vjust = 0)
      )
    print(p1)
  })
  
  
}



for(i in 1){
  agg_png(file=paste0(wdShinyPlots, "plotsSub_CCREM_cohort_DemImportant_onlyInterviewed.png"), width = 3360, height = 2240, res = 300)
  
  grid.arrange(arrangeGrob(pg_list_DemImportant[[i]], pg_list_DemImportant[[i+1]], pg_list_DemImportant[[i+2]], pg_list_DemImportant[[i+3]], pg_list_DemImportant[[i+4]], pg_list_DemImportant[[i+5]], pg_list_DemImportant[[i+6]], pg_list_DemImportant[[i+7]], pg_list_DemImportant[[i+8]], pg_list_DemImportant[[i+9]], pg_list_DemImportant[[i+10]], pg_list_DemImportant[[i+11]], pg_list_DemImportant[[i+12]], pg_list_DemImportant[[i+13]], pg_list_DemImportant[[i+14]], pg_list_DemImportant[[i+15]], nrow = 4))
  dev.off()
}


for(i in seq(1,length(pg_list_Essentials), by = 15)){
  agg_png(file=paste0(wdShinyPlots, "plotsSub_CCREM_cohort_Essentials_",i,"_onlyInterviewed.png"), width = 3360, height = 2240, res = 300)
  
  grid.arrange(arrangeGrob(pg_list_Essentials[[i]], pg_list_Essentials[[i+1]], pg_list_Essentials[[i+2]], pg_list_Essentials[[i+3]], pg_list_Essentials[[i+4]], pg_list_Essentials[[i+5]], pg_list_Essentials[[i+6]], pg_list_Essentials[[i+7]], pg_list_Essentials[[i+8]], pg_list_Essentials[[i+9]], pg_list_Essentials[[i+10]], pg_list_Essentials[[i+11]], pg_list_Essentials[[i+12]], pg_list_Essentials[[i+13]], pg_list_Essentials[[i+14]], nrow = 4))
  dev.off()
}





plc_list_DemImportant <- list()
plc_list_Essentials <- list()



for(i in 1:length(predList_DemImportant)){
  
  message(i)
  
  plc_list_DemImportant[[i]] <- local({
    i <- i 
    pp.data <- predAgeGrp_DemImportant %>% filter(dv_name == name_flag_DemImportant[i], country == country_flag_DemImportant[i])
    p1 <- ggplot() + 
      geom_errorbar(data = pp.data, aes(x=age_grp, ymin = LL, ymax = UL), alpha = 0.2, colour = cols[2]) +
      geom_point(data = pp.data, aes(x=age_grp, y=fit), alpha = 0.1, size = 1, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=age_grp,y = fit), se = F, alpha = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=age_grp,y = LL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=age_grp,y = UL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) +
      ggtitle(paste(pp.data$country[1])) + 
      scale_y_continuous(name="", breaks=seq(-1, 1, 0.4), limits = c(-1,1)) + 
      scale_x_continuous(name = " ",breaks = seq(1,8,1) ,labels = c("15-19", "", "30-39", "", "50-59", "", "70-79", "")) + 
      theme(
        axis.line = element_line(colour = "black"), 
        text = element_text(size=7), 
        axis.text = element_text(size = 7),
        plot.title = element_text(margin = margin(0,0,0,0), size = 7, vjust = -4),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt")
      ) + 
      theme_light() + 
      theme(
        plot.title = element_text(margin = margin(0,0,0,0)),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt"),
        axis.text.x.bottom = element_text(angle = -90, hjust = 0, vjust = 0)
      )
    print(p1)
  })
  
  
}



for(i in 1:length(predList_Essentials)){
  
  message(i)
  
  plc_list_Essentials[[i]] <- local({
    i <- i 
    pp.data <- predAgeGrp_Essentials %>% filter(dv_name == name_flag_Essentials[i], country == country_flag_Essentials[i])
    p1 <- ggplot() + 
      geom_errorbar(data = pp.data, aes(x=age_grp, ymin = LL, ymax = UL), alpha = 0.2, colour = cols[2]) +
      geom_point(data = pp.data, aes(x=age_grp, y=fit), alpha = 0.1, size = 1, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=age_grp,y = fit), se = F, alpha = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=age_grp,y = LL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=age_grp,y = UL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) +
      ggtitle(paste(pp.data$country[1])) + 
      scale_y_continuous(name="", breaks=seq(-1, 1, 0.4), limits = c(-1,1)) + 
      scale_x_continuous(name = " ",breaks = seq(1,8,1) ,labels = c("15-19", "", "30-39", "", "50-59", "", "70-79", "")) + 
      theme(
        axis.line = element_line(colour = "black"), 
        text = element_text(size=7), 
        axis.text = element_text(size = 7),
        plot.title = element_text(margin = margin(0,0,0,0), size = 7, vjust = -4),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt")
      ) + 
      theme_light() + 
      theme(
        plot.title = element_text(margin = margin(0,0,0,0)),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt"),
        axis.text.x.bottom = element_text(angle = -90, hjust = 0, vjust = 0)
      )
    print(p1)
  })
  
  
}



for(i in 1){
  agg_png(file=paste0(wdShinyPlots, "plotsSub_CCREM_lifecycle_DemImportant_onlyInterviewed.png"), width = 3360, height = 2240, res = 300)
  
  grid.arrange(arrangeGrob(plc_list_DemImportant[[i]], plc_list_DemImportant[[i+1]], plc_list_DemImportant[[i+2]], plc_list_DemImportant[[i+3]], plc_list_DemImportant[[i+4]], plc_list_DemImportant[[i+5]], plc_list_DemImportant[[i+6]], plc_list_DemImportant[[i+7]], plc_list_DemImportant[[i+8]], plc_list_DemImportant[[i+9]], plc_list_DemImportant[[i+10]], plc_list_DemImportant[[i+11]], plc_list_DemImportant[[i+12]], plc_list_DemImportant[[i+13]], plc_list_DemImportant[[i+14]], plc_list_DemImportant[[i+15]], nrow = 4))
  dev.off()
}


for(i in seq(1,length(plc_list_Essentials), by = 15)){
  agg_png(file=paste0(wdShinyPlots, "plotsSub_CCREM_lifecycle_Essentials_",i,"_onlyInterviewed.png"), width = 3360, height = 2240, res = 300)
  
  grid.arrange(arrangeGrob(plc_list_Essentials[[i]], plc_list_Essentials[[i+1]], plc_list_Essentials[[i+2]], plc_list_Essentials[[i+3]], plc_list_Essentials[[i+4]], plc_list_Essentials[[i+5]], plc_list_Essentials[[i+6]], plc_list_Essentials[[i+7]], plc_list_Essentials[[i+8]], plc_list_Essentials[[i+9]], plc_list_Essentials[[i+10]], plc_list_Essentials[[i+11]], plc_list_Essentials[[i+12]], plc_list_Essentials[[i+13]], plc_list_Essentials[[i+14]], nrow = 4))
  dev.off()
}



